{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# StackOverflow Problems"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Real-world problems to test your skills on!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from skimage import (filters, io, color, exposure, feature,\n",
    "                     segmentation, morphology, img_as_float)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Parameters of a pill\n",
    "\n",
    "(Based on StackOverflow http://stackoverflow.com/questions/28281742/fitting-a-circle-to-a-binary-image)\n",
    "\n",
    "<img src=\"../../images/round_pill.jpg\" width=\"200px\" style=\"float: left; padding-right: 1em;\"/>\n",
    "Consider a pill from the [NLM Pill Image Recognition Pilot](http://pir.nlm.nih.gov/pilot/instructions.html) (``../../images/round_pill.jpg``).  Fit a circle to the pill outline and compute its area.\n",
    "\n",
    "<div style=\"clear: both;\"></div>\n",
    "\n",
    "*Hints:*\n",
    "\n",
    "1. Equalize (``exposure.equalize_*``)\n",
    "2. Detect edges (``filter.canny`` or ``feature.canny``--depending on your version)\n",
    "3. Fit the ``CircleModel`` using ``measure.ransac``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "image = io.imread(\"../../images/round_pill.jpg\")\n",
    "image_equalized = exposure.equalize_adapthist(image)\n",
    "edges = feature.canny(color.rgb2gray(image_equalized))\n",
    "\n",
    "f, (ax0, ax1, ax2) = plt.subplots(1, 3, figsize=(15, 8))\n",
    "ax0.imshow(image)\n",
    "ax1.imshow(image_equalized)\n",
    "ax2.imshow(edges, cmap='gray');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage import measure\n",
    "from matplotlib.patches import Circle\n",
    "\n",
    "coords = np.column_stack(np.nonzero(edges))\n",
    "\n",
    "model, inliers = measure.ransac(coords, measure.CircleModel,\n",
    "                                min_samples=3, residual_threshold=1,\n",
    "                                max_trials=500)\n",
    "\n",
    "print('Circle parameters:', model.params)\n",
    "\n",
    "row, col, radius = model.params\n",
    "\n",
    "f, ax = plt.subplots()\n",
    "ax.imshow(image, cmap='gray');\n",
    "circle = Circle((col, row), radius=radius, edgecolor='C9', linewidth=2, fill=False)\n",
    "ax.add_artist(circle);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Alternative: morphological snakes\n",
    "\n",
    "**NOTE**: this is expensive to compute, so may take a while to execute"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initial level set\n",
    "pill = color.rgb2gray(image)\n",
    "pill = restoration.denoise_nl_means(pill, multichannel=False)\n",
    "\n",
    "level_set = segmentation.circle_level_set(pill.shape, radius=200)\n",
    "ls = segmentation.morphological_chan_vese(pill, 80, init_level_set=level_set, smoothing=3)\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, figsize=(8, 8))\n",
    "\n",
    "ax.imshow(pill, cmap=\"gray\")\n",
    "ax.set_axis_off()\n",
    "ax.contour(ls, [0.5], colors='r');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Counting coins\n",
    "\n",
    "Based on StackOverflow http://stackoverflow.com/questions/28242274/count-number-of-objects-using-watershed-algorithm-scikit-image\n",
    "\n",
    "Consider the coins image from the scikit-image example dataset, shown below.\n",
    "Write a function to count the number of coins.\n",
    "\n",
    "The procedure outlined here is a bit simpler than in the notebook lecture (and works just fine!)\n",
    "\n",
    "<div style=\"clear: both;\"></div>\n",
    "\n",
    "*Hint:*\n",
    "\n",
    "1. Equalize\n",
    "2. Threshold (``filters.threshold_otsu``)\n",
    "3. Remove objects touching boundary (``segmentation.clear_border``)\n",
    "4. Apply morphological closing (``morphology.closing``)\n",
    "5. Remove small objects (``measure.regionprops``)\n",
    "6. Visualize (potentially using ``color.label2rgb``)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage import data\n",
    "fig, ax = plt.subplots()\n",
    "ax.imshow(data.coins(), cmap='gray');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage import segmentation\n",
    "\n",
    "image = data.coins()\n",
    "equalized = exposure.equalize_adapthist(image)\n",
    "binary0 = equalized > filters.threshold_otsu(equalized)\n",
    "binary1 = segmentation.clear_border(binary0)\n",
    "binary2 = morphology.closing(binary1, morphology.square(3))\n",
    "\n",
    "f, (ax0, ax1) = plt.subplots(1, 2)\n",
    "ax0.imshow(image, cmap='gray')\n",
    "ax1.imshow(edges, cmap='gray');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "labels = ndi.label(binary2)[0]\n",
    "labels_big = morphology.remove_small_objects(labels)\n",
    "print(\"Number of coins:\", len(np.unique(labels_big)[1:]))\n",
    "\n",
    "out = color.label2rgb(labels_big, image, bg_label=0)\n",
    "fig, ax = plt.subplots()\n",
    "ax.imshow(out);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Snakes\n",
    "\n",
    "Based on https://stackoverflow.com/q/8686926/214686\n",
    "\n",
    "<img src=\"../../images/snakes.png\" width=\"200px\" style=\"float: left; padding-right: 1em;\"/>\n",
    "\n",
    "Consider the zig-zaggy snakes on the left (``../../images/snakes.png``).<br/>Write some code to find the begin- and end-points of each.\n",
    "\n",
    "<div style=\"clear: both;\"></div>\n",
    "\n",
    "*Hints:*\n",
    "\n",
    "1. Threshold the image to turn it into \"black and white\"\n",
    "2. Not all lines are a single pixel thick.  Use skeletonization to thin them out (``morphology.skeletonize``)\n",
    "3. Locate all snake endpoints (I used a combination of ``scipy.signal.convolve2d`` [find all points with only one neighbor], and ``np.logical_and`` [which of those points lie on the snake?] to do that, but there are many other ways)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import ndimage as ndi\n",
    "\n",
    "image = color.rgb2gray(io.imread('../../images/snakes.png'))\n",
    "\n",
    "# Reduce all lines to one pixel thickness\n",
    "snakes = morphology.skeletonize(image < 1).astype(np.uint8)\n",
    "\n",
    "# Find pixels with only one neighbor\n",
    "neighbor_kernel = np.array([[1, 1, 1],\n",
    "                            [1, 0, 1],\n",
    "                            [1, 1, 1]])\n",
    "num_neighbors = ndi.convolve(snakes, neighbor_kernel,\n",
    "                             mode='constant')\n",
    "corners = (num_neighbors == 1) & snakes\n",
    "\n",
    "# Those are the start and end positions of the segments\n",
    "rr, cc = np.nonzero(corners)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.imshow(image, cmap='gray')\n",
    "ax.scatter(cc, rr)\n",
    "ax.set_axis_off()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# M&Ms\n",
    "\n",
    "How many blue M&Ms are there in this image (`../../images/mm.jpg`)?\n",
    "\n",
    "<img src=\"../../images/mm.jpg\" width=\"400px\"/>\n",
    "\n",
    "Steps:\n",
    "\n",
    "1. Denoise the image (using, e.g., `restoration.denoise_nl_means`)\n",
    "2. Calculate how far each pixel is away from pure blue\n",
    "3. Segment this distance map to give a \"pill mask\"\n",
    "4. Fill in any holes in that mask, using `scipy.ndimage.binary_fill_holes`\n",
    "5. Use watershed segmentation to split apart any M&Ms that were joined, as described in http://scikit-image.org/docs/dev/auto_examples/segmentation/plot_watershed.html\n",
    "\n",
    "*Alternative approach:*\n",
    "\n",
    "- http://scikit-image.org/docs/dev/user_guide/tutorial_segmentation.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "image = img_as_float(io.imread('../../images/mm.jpg'))\n",
    "blurred = restoration.denoise_nl_means(image, multichannel=True)\n",
    "blue_pill_color = [0, 0, 1]  # rough approximation; good enough\n",
    "\n",
    "blue_dist = np.linalg.norm(blurred - [0, 0, 1], axis=2)\n",
    "blue_mask = blue_dist < np.percentile(blue_dist, 5)\n",
    "blue_mask = ndi.binary_fill_holes(blue_mask)\n",
    "\n",
    "plt.imshow(blue_mask, cmap='gray');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "distance = ndi.distance_transform_edt(blue_mask)\n",
    "local_maxi = feature.peak_local_max(\n",
    "    distance, indices=False, footprint=np.ones((5, 5)), labels=blue_mask\n",
    ")\n",
    "markers = ndi.label(local_maxi)[0]\n",
    "labels = segmentation.watershed(-distance, markers, mask=blue_mask)\n",
    "\n",
    "plt.imshow(labels, cmap='tab20');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"There are {} blue M&M's\".format(np.max(labels)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Viscous fingers\n",
    "\n",
    "Based on StackOverflow: http://stackoverflow.com/questions/23121416/long-boundary-detection-in-a-noisy-image\n",
    "\n",
    "<img src=\"../../images/fingers.png\" width=\"200px\" style=\"float: left; padding-right: 1em;\"/>\n",
    "\n",
    "Consider the fluid experiment on the right.  Determine any kind of meaningful boundary in this noisy image.\n",
    "\n",
    "<div style=\"clear: both;\"></div>\n",
    "\n",
    "*Hints:*\n",
    "\n",
    "1. Convert to grayscale\n",
    "2. Try edge detection (``feature.canny``)\n",
    "3. If edge detection fails, denoising is needed (try ``restoration.denoise_tv_bregman``)\n",
    "4. Try edge detection (``feature.canny``)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage import restoration, color, io, feature, morphology\n",
    "\n",
    "image = color.rgb2gray(img_as_float(io.imread('../../images/fingers.png')))\n",
    "denoised = restoration.denoise_nl_means(image, h=0.06, multichannel=False)\n",
    "edges = feature.canny(denoised, low_threshold=0.001, high_threshold=0.75, sigma=1)\n",
    "\n",
    "fig, (ax0, ax1) = plt.subplots(1, 2, figsize=(20, 15))\n",
    "ax0.imshow(denoised, cmap='gray')\n",
    "ax1.imshow(edges, cmap='gray', interpolation='lanczos')\n",
    "for ax in (ax0, ax1):\n",
    "    ax.set_axis_off()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.6rc1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
